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OBJECTIVE 


Develop  a mathematical  model  for  the  wideband  HF  channel  and 
the  software  implementation  of  this  model.  The  model  must  be 
valid  for  carrier  frequencies  in  the  2-32  MHz  range,  for  band- 
widths  up  to  100  kHz,  and  for  ground  ranges  up  to  1000  km. 

RESULTS 

A simplified  model  of  the  ionosphere  is  employed  to  deter- 
mine the  various  skywave  returns  and  the  effect  these  returns  have 
on  the  overall  received  signal.  The  groundwave  is  modeled  as  a 
delayed,  attenuated  version  of  the  transmitted  signal.  Inputs  to 
the  simulator  consist  of  those  basic  physical  quantities  that 
characterize  the  ionosphere  and  those  that  characterize  the 
transmitter-receiver  environment;  values  of  delay,  attenuation, 
doppler  shift  and  spread,  and  dispersion  are  calculated  internally 
for  each  component  of  the  received  signal.  A FORTRAN  version  of 
the  simulator  has  been  written  and  tested  for  several  typical 
channels . 

RECOMMENDATIONS 

Make  further  refinements  in  the  mathematical  model  in  order 
to  expand  the  scope  and  accuracy  of  the  simulator.  Of  particular 
value  in  this  regard  would  be  the  inclusion  of  a deviative  absorp- 
tion term  for  the  skywave  returns  and  the  inclusion  of  an  attenua- 
tion term  for  propagation  of  the  groundwave  over  a rough,  lossy 
ocean  surface. 

Investigate  the  relationship  between  the  simulator  inputs 
and  the  actual  observed  returns  as  measured,  for  example,  by  an 
ionosonde.  A knowledge  of  this  relationship  would  provide  further 
verification  of  the  mathematical  model  as  well  as  a set  of  inputs 
that  correspond  to  known  operating  conditions. 

Finally,  in  view  of  the  meager  data  base  available  in  the 
wideband  case,  any  measurement  program  aimed  at  enlarging  this 
data  base  would  represent  a significant  contribution  to  the  study 
and  modeling  of  HF  channels. 
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GLOSSARY 


English  Letter  Symbols 


a 

Multiplicative  constant  used  to 
generate  random  sequence  {Y^^ 

(n) 

Impulse  response  of  filter  used  to 
ate  approximate  random  gain 

A 

Multiplicative  constant  for  doppler 
shift 

\ 

Gain  of  groundwave 

A. 

D 

RMS  gain  of  j-th  ionosphere  path 

A(x) 

= 2 G tan  6 

B 

Multiplicative  constant  for  doppler 
spread 

B 

Bearing  of  receiver  as  seen  from 

10 

transmitter  in  magnetic  coordinates 

BW 

Two-sided  bandv/idth  of  transmitted 
signal 

C 

Velocity  of  light 

C_  (T) 

Correlation  function  of  G^  (t) 

J 

(T) 

Correlation  function  of  &.(t) 

D 

D 

Correlation  function  of  (t)  and 

D 

virtual  path  length 

D 

g 

Delay  of  groundwave 

= f . 

1 

c 3 

3 

f 

Frequency 

f 

c 

Carrier  frequency 

VI 


GLOSSARY  (Continued) 


f (x) 

^pE 


(o 

pF 


f (x) 


G = G(f) 


Doppler  reference  frequency 
Gyrof requency  of  d region 
Gyrof requency  of  E region 

Gyrof requency  of  F region 

General  expression  for  penetration  fre 
quency  of  E region  (ordinary  or  extra- 
ordinary wave) 

General  expression  for  penetration  fre 
quency  of  F region  (ordinary  or  extra- 
ordinary wave) 

Penetration  frequency  of  E region  - 
ordinary  wave 

Penetration  frequency  of  E region  - 
extraordinary  wave 

Penetration  frequency  of  F region  - 
ordinary  wave 

Penetration  frequency  of  F region  - 
extraordinary  wave 

Group  height  of  signal  path,  oblique 
incidence 


Gj  (t) 

Sj  (t) 

G^(6) 

G,i,(e) 


Random  gain  of  j-th  path 

Approximation  for  random  gain  of  j-th 
path 

Gain  of  receiving  antenna 

Gain  of  transmitting  antenna 

Group  height  of  signal  path,  vertical 
incidence 


h 


Height  above  earth 


Height  of  peak  electron  density  of 
D region 


^E 


Height  of  peak  electron  density  of 
E region 
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GLOSSARY  (Continued) 


Height  of  peak  electron  density  of 
F region 


h 


r 


Height  at  which  refractive  index 
equals  zero  at  vertical  incidence 


h.  (t) 


Distortion  impulse  response  of  j-th 
path 


h.  (t) 


2H. 


Modified  distortion  impulse  response 

=T  h. (t) 
o : 

Semi thickness  of  E region  electron 
density 


Hj  (f) 
I 

I(x) 

I (x) 
P 


Semi thickness  of  F region  electron 
density 

Fourier  transform  of  h^ (t) 

Angle  of  inclination  of  earth's 
magnetic  field 

Cannonical  expression  representing 
general  impulse  response  hj (t) 

Integrals  used  in  series  representa- 
tion of  I (x) 


k 

K 


^bs 

^ant 

Ld 

M 

N=N(h) 


Sunspot  number  multiplier  for  non- 
deviative  absorption  loss 

Multiplicative  constant  for  non- 
deviative  absorption  loss 

Constant  used  in  the  generation  of 

^j(nT^) 

Absorption  loss 

Antenna  pattern  loss 

Spreading  loss 

Total  loss  of  j-th  path 

Length  of  pseudo  random  sequence 

Electron  density 

Maximum  electron  density  of  E region 
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GLOSSARY  (Continued) 


N. 


N 


N . 
D 


Pp  (f) 


r(t) 

R 


s (t) 


Sj  (t) 


§.  (t) 


S 

S(f) 

t 


Ui,U2 


(t)  (t) 


Wj  (1)  ,Wj  (2)  ,Wj  (3) 


^^n> 


Maximum  electron  density  of  F region 


Integer  multiple  of  T nearest  the 
groundwave  delay  ° 


Integer  multiple  of  T nearest  the 
signal  delay  ° 


Power  spectral  density  of  (t) 


Power  spectral  density  of  (t) 


(1) 


Power  spectral  density  of  U ' (t)  and 


(t) 


Received,  complex,  baseband  signal 


Distance  between  transmitter  and 
receiver 


Transmitted,  complex,  baseband  signal 


Distorted,  complex  signal  returned 
from  ionosphere  via  j-th  path 


Distorted,  complex  signal  returned  from 
ionosphere  via  j-th  path  after  removal 
of  phase  and  delay  components: 


s.  (t)  = exp  [-2iTiD  . ] s.(t-D!^S 

3 ^ 3 3 3 


Sunspot  number 

Fourier  transform  of  s(t) 

Time 


Sampling  interval  (one  half  the 
Nyquist  interval  of  the  transmitted 
signal) 


Uniform  [0,1]  random  variables 


Real,  independent,  zero  mean  identically- 
distributed,  Gaussian  process  used  to 
define  G . (t) 


Independent,  complex,  Gaussian  random 
variables  used  in  generation  of  G^ (nT^) 


= cos  9 

Sequence  of  pseudo  random  integers 
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Greek  Letter  Symbols 


a 

e 

Y 


Yj(l),  Yj(2) 


(2) 

(3) 


A . 
3 


(2) 


A . 
3 


(3) 


Doppler  shift  exponent 

Doppler  spread  exponent 

Exponent  of  solar  zenith  angle  term  in 
expression  for  nondeviative  absorption 

Constants  used  in  the  generation  of 

a . (nT„) 

General  amplitude  distortion  coefficient 

General  phase  distortion  coefficient 

Normalized  amplitude  distortion 
coefficient  of  j-th  path 

Normalized  phase  distortion  coefficient 
of  j-th  path 

Normalized  frequency  variable  = f/BW 


0 


K 

X 


Ray  angle  of  signal  path  measured  from 
vertical 

= (e^/4TT^eQm)  = 80.5 
Geographic  longitude 


X 

o 


Geographic  longitude  of  north  magnetic 
pole 


X . (1) , X . (2)  Constants  used  in  the  generation  of 

3 3 

G . (nT^) 

U Refractive  index 


V 


General  doppler  shift 


V . 

3 


Doppler  shift  of  j-th  path 
Reference  doppler  shift 


i 


X 


Cj  (n) 


2o 


2a  . 
3 

2o_ 


Tj  (f) 


T . 

D 


(k) 


<J> 

X 


GLOSSARY  (Continued) 


Complex  Gaussian  random  variable 

Independent,  identically-distributed 
Gaussian  sequence 

Real  part  of  (n) 

Imaginary  part  of  (n) 

Parameter  used  in  fitting  Cp<  (t)  to 
C_  (T)  ^3 

General  doppler  spread 

Doppler  spread  of  j-th  path 

Reference  doppler  spread 

General  path  delay 

Delay  of  3-th  path  as  a function  of 
frequency 

k-th  derivative  of  delay  of  j-th  path 
evaluated  at  the  carrier  frequency  f^ 

Geographic  latitude 

Geographic  latitude  of  north  magnetic 
pole 

Magnetic  latitude 
Solar  zenith  angle 

Angle  between  earth's  magnetic  field  and 
wave  incident  upon  ionosphere 

Angle  between  earth's  magnetic  field 
and  wave  reflected  from  ionosphere 
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INTRODUCTION 


The  use  of  channel  simulators  to  study  the  performance  of 
communication  systems  and  to  assist  in  their  design  has  become 
increasingly  popular  in  recent  years.  By  employing  a simulator, 
it  is  possible  to  reproduce  given  channel  characteristics  as 
desired  and  to  compare  different  system  designs  in  identical 
environments.  Moreover,  the  cost  of  simulation  is  usually  much 
less  than  the  cost  of  "on  the  air"  measurements  and  it  is  quite 
reliable.  A variety  of  HF  channel  simulators  have  been  pro- 
posed and  studied  in  the  past  [1-4]*.  Sailors  and  Hill  [5]  have 
written  an  excellent  survey  of  the  field  and  have  made  several 
recommendations  for  future  work.  Some  of  these  recommendations  are 
embodied  in  the  present  report. 

In  the  following  the  development  of  a simulator  for  the  wide- 
band HF  channel  is  described.  This  simulator  is  designed  for 
carrier  frequencies  in  the  2-32  MKz  range  and  it  can  accommodate 
bandwidths  up  to  100  kHz  and  ground  ranges  up  to  1000  km.  The 
simulator  models  the  groundwave  and  all  ionosphere  returns  (both 
single  and  multiple  hops) , but  does  not  include  additive  noise  or 
signal  distortion  produced  in  the  transmitter  or  receiver.  It  takes 
into  account  all  ionospheric  phenomena  that  have  a significant 
effect  on  wideband  communications.  Many  of  these  phenomena,  such 
as  delay  differences  between  the  ordinary  and  extraordinary  waves 
and  signal  distortion  caused  by  variations  of  delay  with  frequency, 
have  been  neglected  in  the  earlier  simulators,  resulting  in  severe 
bandwidth  limitations.  The  inputs  required  are  the  basic  para- 
meters that  characterize  the  physical  ionosphere  (height,  thickness, 
and  maximum  electron  density  of  the  E and  F regions,  sunspot 
number,  time  of  day,  etc.)  and  those  that  characterize  the  trans- 
mitter and  receiver  (carrier  frequency,  range,  geographic  location, 
antenna  patterns,  etc.).  Since  these  input  parameters  are  the  ones 
that  are  actually  known  in  practice  it  is  a simple  matter  to  vary 
them  and  observe  the  effect  they  have  on  system  performance.  This 

*References  are  listed  in  Section  VI. 
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capability  does  not  exist  in  the  previous  simulators  which  require 
as  inputs  not  the  basic  quantities  mentioned  above  but  those 
channel  parameters  such  as  delay,  attenuation,  doppler  shift  and 
doppler  spread  derived  from  them. 

In  the  course  of  developing  the  simulator  a number  of  simpli- 
fications and  approximations  have  been  made  to  obtain  tractable 
results.  The  most  important  of  these  assumptions  concerns  the 
electron  density  profile  and  the  absorption  mechanism.  It  is 
recognized  that  the  ionosphere  model  used  does  not  provide  a detailed 
physical  characterization  of  the  ionosphere;  however,  it  does  pro- 
vide an  "order  of  magnitude"  characterization  suitable  for  communi- 
cation purposes.  The  simulator  is  divided  into  two  parts.  The 
first  part  is  concerned  with  the  mathematical  model  of  the  iono- 
sphere and  the  computation  of  the  channel  parameters  from  the  basic 
inputs.  The  second  part  uses  these  parameters  to  determine  the 
channel  response  for  arbitrary  inputs.  Both  parts  of  the  simulator 
are  described  in  detail  in  the  following  sections. 
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II.  THE  HF  CHANNEL 


The  signal  received  via  the  HF  channels  considered  in  this 
report  is  the  sum  of  the  groundwave  and  the  various  ionosphere 
returns.  Throughout  we  shall  assume  that  the  groundwave  is  an 
attenuated,  delayed  version  of  the  transmitted  signal  but  is 
otherwise  undistorted.  The  ionosphere  returns  also  suffer  attenu- 
ation and  delay  but  in  addition  are  subject  to  fading  and  distor- 
tion. There  is  also  time  dispersion  (multipath)  produced  by  the 
different  delays  associated  with  each  ionosphere  return.  The 
purpose  of  this  section  is  to  present  some  definitions  that  are 
used  in  the  sequel  and  to  introduce  the  mathematical  model  for 
the  HF  channel. 

The  complex  transmitted  baseband  signal  will  be  denoted  by 
s(t)  and  its  Fourier  transform  by  S(f).  Thus 

00 

s{t)  =J'  S(f)  exp  [2Trift]df.  (1) 

— oo 

Distorted  versions  of  this  signal  will  be  returned  from  the  iono- 
sphere for  both  the  E region  and  the  F region  (low  ray  and  high 
ray) . There  will  be  a set  of  these  returns  for  the  ordinary  and 
extraordinary  magnetoionic  components  and  for  as  many  multiple-hop 
paths  as  are  present.  Although  many  returns  are  possible,  all  but 
a few  of  these  experience  a large  attenuation  and  the  number 
of  effective  returns  is  usually  small.  In  the  following  we  will 
denote  the  number  of  effective  ionosphere  returns  by  N.  This 
number  is  determined  by  excluding  all  returns  that  suffer  an 
attenuation  greater  than  some  predetermined  value. 

Distortion  of  signals  reflected  from  the  ionosphere  has  been 
studied  by  several  authors,  who  have  attempted  to  characterize 
the  impulse  response  for  a simplified  model  of  the  ionosphere  [6-10] . 
Unfortunately,  these  results  have  been  derived  using  ionosphere 
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models  that  are  too  idealized  for  our  purposes'*'  and  the  impulse 

2 

responses  obtained  are  too  complicated  for  easy  implementation  . 
In  order  to  obtain  tractable  expressions  for  the  distortion  v/e 
shall  assume  that  it  is  caused  primarily  by  the  fact  that  the 
delay  is  a function  of  frequency  and  is  not  constant  across  the 
signal  bandwidth.  This  characterization  of  the  distortion  mech- 
anism allows  us  to  make  a number  of  approximations  which  result 
in  relatively  simple  expressions  for  the  distorted  signal.  For 
obvious  reasons  we  shall  refer  to  this  particular  distortion  as 
delay  distortion. 

Let  f^  denote  the  carrier  frequency  and  let  t . (f ) denote  the 

3 

delay  of  the  j-th  return  as  a function  of  frequency.  It  follows 
from  (1)  that  the  distorted  baseband  signal  returned  from  the 
ionosphere  is  given  by 


s . (t)  =/  S(f)  exp  [-2T7i  (f  +f ) -T  . (f  +f  )]*exp  [2TTift]  df. 

3 •/  3 

— CO 

j — 1 , • • • , N 


(2) 


Vie  will  characterize  this  distortion  by  the  coefficients  in  the 
power  series  expansion  of  the  delay  about  the  carrier  frequency. 
This  series  appears  as 


Efk  ^ (k) 

--  "j 

^ k=0  k! 


j— 1 / • • • /N 


:3) 


^They  allow  for  reflection  only  from  a single  ionospheric  region 
and  consider  such  simple  electron  density  profiles  as  a constant 
and  a linear  function  of  height. 

2 

In  one  case  the  impulse  response  is  given  as  the  solution  of  an 
integral  equation. 


4 


with 


D ^ c 


Substituting  this  series  into  the  exponential  delay  term  in  (2)  we 
obtain 


exp  [-2Tri  (f^+f)>T  j (f^+f )]  = exp 


where  we  have  defined 


r-2TTi  ^ Dj 

L k=0 


D.<°)  = f T.<°) 
D CD 


j=l, • • • ,N 


j =1 , • • • N 


*-  -•  j=l,--.M. 


Thus,  for  Sj (t)  we  have 


UU  oo 

5j  (t)  = exp^-2TriDj  j J'  S (f)exp  p2TTi  f^D^ 

-00  L 1,=?  J 


■exp  ^2iTif  (t-Dj  )j  df . 


j=l, • • • ,N 


For  computational  reasons  it  is  convenient  to  separate  the 
phase  and  delay  terms  from  the  actual  distortion  of  the  signal. 
To  this  end  we  define  a function  s^ (t)  as  follows 


s.(t)  = exp  [-2TTiD  . s.(t-D.^^^) 

D D D D 


j=l , • * • N 
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where 


5.(t)  = f S (f ) H . (f ) exp  [ 2TTif  t]  dt 
3 3 


(10) 


j=l, • • • ,N 


with 


oo 

Hj(f)  = exp|-2iri  ^ , j = l,'*-N. 

^ k=2  -* 


(11) 


In  practice  it  is  convenient  to  determine  (t)  by  performing  a 
time  domain  convolution  rather  than  by  evaluating  the  Fourier 
transform  (10) . If  the  impulse  response  determined  by  (f)  is 
denoted  by  h^ (t) , we  have 

OO 

h.(t)  = / H.(f)  exp  [27Tift]df,j=l,  • • - N.  , 

D D 

— 00 

and  the  desired  convolution  appears  as 


(12) 


Sj  (t) 


= / 


s (t-u) hj 


(u)du,  j=l. 


•N, 


(13) 


A discussion  of  the  approximations  used  in  the  numerical  evaluation 

/\ 

of  hj (t)  and  s^ (t)  is  presented  in  Section  IV. 

In  addition  to  the  above  distortion,  the  ionosphere  returns 
are  subject  to  attenuation  and  random  fading.  These  two  phenomena 
are  included  by  multiplying  the  delayed,  distorted  signal  by  a 
random  gain,  which  we  denote  by  (t) . We  will  assume  that  (t) 
is  a complex,  zero-mean,  Gaussian  process  with  identically-distributed 
real  and  imaginary  parts.  This  is  the  Rayleigh  fading  model  [11] . 
While  the  assumption  of  Rician  fading  [12]  would  provide  somewhat 
more  flexibility  than  the  Rayleigh  assumption,  the  available 
evidence  [12-14]  indicates  that  the  Rayleigh  model  is  adequate 
for  our  purposes.  The  random  gains  associated  with  different  returns 
are  assumed  to  be  independent. 


ii 
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zero  mean,  identically-distributed,  stationary  Gaussian  processes 
with  correlation  function 

Cy  (t)  = E j (t+T)U^^^ j (t)  = E U j {t+T)U j (t) . (14) 

j=l, • • -N, 

The  corresponding  power  spectral  density  is  given  by 


Py  (f)  = /*  C (x)  exp[-2TrifT]dT  , j=l,  • • • ,N. 

3 _«  3 

We  define  the  random  gain  in  terms  of  (t)  and  (t)  as 


follows 


Gj(t)  = [Uj^^^(t)  + i Uj^^Nt)]-exp[2TTiVjt] 

j=l, • • -N. 

where  is  the  doppler  shift  for  the  j-th  ionosphere  return. 

The  correlation  and  power  spectral  density  of  G.(t)  are  easily 

( 1 )3  ( 2 ) 

determined  using  the  assumed  properties  of  and  (t)  . 

Thus 


C (t)  = EG  . (t+T)G  . (t)  = 2C  (t)*  exp[2TTiv  .t] 
Gj  3 3 Uj  3 


j=l, • • ’N. 


P_  (f)  = 2P  (f-v.)  , j=l,---N.  (18) 

Gj  Uj  3 

Following  Watterson  [1]  we  assume  that  the  power  spectral 
density  P_  (f)  has  a Gaussian  shape.  While  there  is  no  conclusive 

experimental  evidence  for  this  assumption,  it  is  a convenient  model 
for  spectra  that  are  unimodal  and  that  decrease  rapidly.  Specifically 
we  let 
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(19) 


2 2 2 2 

C (t)  = A.  exp[-2Ti  o.  T + 2TTiv.T] 
Gj  3 3 3 


j=l» 


•N 


Here  is  the  rms  gain  ( 0 £ A^  £l)  of  the  j-th  return  and  2a ^ is 
the  doppler  spread  of  the  return.  The  power  spectral  density 
associated  with  this  correlation  function  is 


3 


r 1 

-< exp i , 3 = 1,  ••• 

2tt  a . I 2a . J 

3 3 


A 

/n 


N. 


(20) 


Finally,  we  must  add  the  groundwave  to  the  ionosphere  returns, 
The  groundwave  is  assumed  to  experience  only  delay  and  to  have  a 
non-random  gain.  If  A denotes  the  gain  (0  £ A £ D and  D 
denotes  the  delay,  the  contribution  to  the  received  signal  from 
the  groundwave  is  given  by  A^  s(t-Dg).  The  gain  is  a complicated 
function  of  range,  sea  state  and  frequency.  Graphs  of  this  func- 
tional relationship  are  given  in  [15] . 

Summarizing  the  above,  the  complex,  baseband  received  signal 
can  be  expressed  as 

N 

\ lr^\  ^ M ^ 

(211 

j J J J 

j=l 


r(t)  = A s(t-D  ) + exp  [-2TTiD  . G . (t)  s . (t-D  . . 

9 9 3 3 3 3 


A block  diagram  illustrating  (21)  is  presented  in  Figure  1. 

In  the  following  sections,  the  quantities  needed  in  the  evalu- 
ation of  r(t)  are  derived  and  the  techniques  employed  in  its  numeri- 
cal evaluation  are  described. 
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III.  MATHEMATICAL  MODEL  OF  THE  IONOSPHERE 


In  this  section  we  introduce  a mathematical  model  of  the 
ionosphere  and  discuss  how  the  channel  parameters  used  in  the 
signal  analysis  part  of  the  simulator  are  determined.  Our  goal 
is  to  obtain  a model  of  the  ionosphere  suitable  for  use  in  mid- 
latitude ocean  environments.  Differences  in  transmission  between 
daytime  and  nighttime  are  included  but  we  do  not  include  those 
effects  related  to  the  rapidly  changing  electron  densities  that 
occur  during  sunrise  and  sunset.  Daytime  splitting  of  the  F region 
into  an  Fj^  and  an  F2  region  is  omitted.  In  addition,  neither 
transmission  via  sporadic  E or  spread  F modes  nor  M-  or  N-type 
reflections  are  included.  In  the  case  of  multiple  hops  we  will 
assume  that  the  ocean  surface  acts  as  a perfect  reflector. 

As  a notational  convenience  in  the  following,  we  will  drop 
the  subscript  j from  all  quantities.  However,  it  should  be  kept 
in  mind  that  each  of  these  quantities  is  associated  with  a par-  [‘ 

ticular  return,  and  it  must  be  identified  with  the  appropriate 
subscript  in  the  signal  analysis  part  of  the  simulator. 

A.  Delay  and  Delay  Derivatives  j 

To  determine  the  group  delay  of  the  signal  reflected  from  ' 

le  ionosphere,  we  make  use  of  the  following  figure. 


Figure  2.  Virtual  and  Actual  Signal  Paths 
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Here  the  dashed  line  is  the  actual  ionospheric  path  of  the  signal 
between  the  transmitter  and  the  receiver  and  the  solid  line  is  the 
virtual  path.  The  ground  range  between  transmitter  and  receiver  is 
denoted  by  R,  the  virtual  path  length  by  D,  the  ray  angle  measured 
from  the  vertical  by  6 and  the  group  height  by  G.  Assuming  a flat 
earth  and  ionosphere,  we  have  by  simple  geometry 

D . — = S (22) 

where  we  have  defined  x = cos  9.  If  c denotes  the  velocity  of 
light,  the  signal  delay  is  given  by  the  virtual  path  length  divided 
by  c 


c(l-x^) 

It  remains  to  determine  the  ray  angle  6. 

For  vertical  incidence  (9=0) , the  group  height  is  given  by 
the  integral  [16] 


G = G (f)  = /*  — 

V V ,T  y 


i 


I 


! 


where  y is  the  refractive  index  and  h^  is  the  height  at  which  y=0. 
In  the  absence  of  both  collisions  and  the  earth's  magnetic  field, 
the  refractive  index  is  given  by  [16] 

y^=  (1-K  45-)  (25) 

f 

where  N = N(h)  is  the  electron  density  in  m ^ , f is  the  frequency 
2 2 

in  Hz  and  < = (e  /4it  t^m)  80.5.  The  group  height  for  oblique 
incidence  is  the  same  as  for  vertical  incidence  but  we  must  replace 
f by  f cos  9 [16],  That  is  G = G(f)  = G^(f  cos  9).  Once  the 
group  height  has  been  determined,  we  appeal  to  Figure  2 to  obtain 
the  relationship 


Defining 


we  see  that,  for  a specified  range  R,  we  must  solve  the  trans- 
cendental equation 

R = A(x)  (28) 


for  X and  then  substitute  this  into  (23)  to  obtain  the  delay. 

We  seek  an  approximation  to  the  electron  density  for  which  the 
group  height,  G,  can  be  expressed  in  a simple  closed  form.  Our 
main  concern  here  is  v/ith  reflections  from  the  E and  F regions. 

At  nighttime  we  assume  that  the  electron  density  profile  can 
be  represented  by  two  distinct  parabolas  as  shown  in  Figure  3. 

-3 

Let  N_  be  the  peak  electron  density  of  the  E region  in  m , let 

hi 

hg  be  the  height  of  the  peak  electron  density  of  the  E region  in 
km,  and  let  2H_  be  the  semithickness  of  the  E region  in  km; 

hi 

similar  definitions  are  made  for  the  F region.  Then  the  night- 
time electron  density  is  given  by 


h-h„  .2 


N(h)  = 


h<h^-2HE 

, hg-2Hg£h£hp+2Hg 
h^+2H^<h<hp-2Hp 
, hp-2Hp<h<hp+2Hp 
, hp+2Hp<h  . 


In  the  daytime  the  electron  density  profile  is  also  assumed 
to  be  parabolic  but  allowance  is  made  for  ionization  between  the 
E and  F regions.  The  electron  density  between  the  two  regions  is 
assumed  to  be  constant  and  equal  to  the  peak  value  of  the  E region. 
This  is  illustrated  in  Figure  4.  Using  the  same  definitions  as 
before,  the  daytime  electron  density  is  given  by 


12 


ELECTRON  DENSITY 


Figure  4.  Daytime  Electron  Density  Pi-otile 
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N(h)  =V 
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r /h-hE  ( 


h-2Hg<h<h^ 


h^<h<h^+2Hj,(l-  (30) 

h ^ 


h +2H^<h. 
F F 


Note  that  the  constants  appearing  in  (24)  and  (30)  will  be  differ- 
ent owing  to  the  different  heights  and  ionization  densities  for 
nighttime  and  daytime  conditions. 

A quantity  of  interest  in  the  present  analysis  is  the  pene- 
tration frequency  of  each  region.  This  is  defined  as  the  plasma 
frequency  at  the  peak  electron  density  [16].  We  have  for  the  E 
and  F regions 


Ola) 


(31b) 


where  »c  = (e  /4tt  ~ 80.5  and  the  two  penetration  frequencies  are 

in  Hz.  The  physical  interpretation  of  fp^  and  f^^  is  that  they 
represent  the  maximum  frequencies  at  which  signals  will  be  reflected 
from  the  respective  regions  at  vertical  incidence. 

Returning  to  the  determination  of  the  group  height,  G,  it  is 
a relatively  simple  matter  to  evaluate  the  integral  (24)  for  either 
of  the  above  electron  densities  [17] . The  details  of  this  evaluation 
are  presented  in  Appendix  A. 

In  the  foregoing  the  effect  of  the  earth's  magnetic  field 
has  been  neglected,  but  we  take  this  to  be  an  adequate  approxima- 
tion for  the  ordinary  wave.  For  the  extraordinary  wave,  however, 
the  penetration  frequencies  are  shifted  relative  to  those  for  the 
ordinary  wave,  and  it  is  essential  to  take  this  into  account. 
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Denoting  the  penetration  frequencies  of  the  ordinary  and  extra- 
ordinary waves  by  the  superscripts  "o"  and  "x"  we  have  [16] 


2 ^HE 


2 ^HF 


f ^ + 4 (f  ^ 
HE  ' pE  ’ 


(32a) 


f„„^  +4  {f^° 
HF  pF 


o)  . 2' 
F ' 


{32b) 


Here  f 


and  f 


^pp  are  defined  by  (31)  , and  f^^  and  f^^p  are  the 

gyrofrequencies  associated  with  the  earth's  magnetic  field  at  the 
levels  of  maximum  ionization  density  in  the  E and  F regions.  For 
the  gyrofrequencies  in  mHz  we  will  use  the  approximate  expressions 
[16] 


^HE  = 


fHF  = 0.87 


' 6370 

6370  + h. 


6370 
6370  + 


^f) 


(1  + 3 sin^  4))^ 


(1  + 3 sin^  <I>) 


(33a) 


(33b) 


where  <I>  is  the  latitude  expressed  in  magnetic  coordinates.  If  (f) 
and  A are  respectively  the  geographic  latitude  and  longitude,  and 
(j)^  and  A^  represent  the  latitude  and  longitude  of  the  north  magne- 


tic pole  (4)  = 78.3°N,  A 

^ o ' o 


69°W) , then  [16] 


sin  4>  = sin  4)  sin  4>  + cos  4>  cos  4>  cos  (A-A  ) . 

T t-o  ^^0  0 


We  assume  that,  to  calculate  delays  for  the  extraordinary  wave, 

(x)  (x) 

we  may  use  the  expressions  in  Appendix  A but  with  fpg  fpp 

substituted  for  the  ordinary  wave  penetration  frequencies.  Other- 
wise the  calculations  of  delay  are  assumed  to  take  the  same  form 
as  those  for  the  ordinary  wave. 
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To  characterize  delay  distortion  we  need  the  derivatives  of 
the  delay  evaluated  at  the  carrier  frequency.  Actually,  only  the 
first  three  derivatives  are  required  for  the  approximation  dis- 
cussed in  Section  IV.  The  determination  of  these  three  deriva- 
tives is  straightforward  although  the  resulting  expressions  are 
quite  complicated.  They  are  listed  in  Appendix  A. 

For  multiple  hop  transmissions  we  assume  that  the  ionosphere 
parameters  are  identical  for  each  hop.  In  the  two-hop  case  the 
foregoing  analysis  is  applied  with  the  original  range  multiplied 
by  1/2.  The  angle  is  then  correct  as  calculated  but  the  path 
length,  delay  and  delay  derivatives  must  be  multiplied  by  2. 

Higher  numbers  of  hops  are  treated  in  a similar  manner  using 
appropriate  multipliers. 

B.  Attenuation 

Three  main  sources  of  attenuation  are  included  in  the  present 

model.  They  are:  (1)  nondeviative  absorption,  (2)  antenna  gain, 

2 

and  (3)  spreading  (1/r  loss) . Other  factors  such  as  focusing 
and  polarization  mismatch  will  also  affect  the  received  signal 
strength,  but  these  are  assumed  to  be  of  secondary  importance  and 
they  have  accordingly  been  omitted. 

B. 1 Absorption 

Absorption  is  assumed  to  take  place  only  in  the  D region  of 
the  ionosphere  and  only  in  the  daytime.  The  absorption  is  assumed 
to  be  nondeviative,  i.e.,  it  takes  place  in  a region  where  the 
refractive  index  is  near  unity.  The  attenuation  will  be  charact- 
erized by  the  expression  [16] 
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(35) 


Here  is  the  attenuation  due  to  absorption  in  decibels  per  hop, 

0 is  the  ray  angle  measured  from  the  vertical,  S is  the  sunspot 
number,  x is  the  solar  zenith  angle,  f is  the  frequency  in  mHz,  f 

HD 

is  the  gyrofrequency  at  the  D region  level  in  mHz  and  [resp.  , 
is  the  angle  between  the  earth's  magnetic  field  and  the  incident 
[resp.,  reflected]  wave.  The  plus  sign  is  used  for  the  ordinary  wave 
and  the  minus  sign  is  used  for  the  extraordinary  wave. 

The  relative  importance  of  the  sunspot  number  and  solar  zenith 
angle  depends  upon  the  constants  K,  k and  y;  following  Davies  [16, 
p.  235]  we  will  set  K = 215,  k = 0.0035  and  y = 0.75.  The  gyro- 
frequency,  is  given  by  [16] 

^HD  " (557^/  I 

where  h^^  is  the  height  of  the  D region  in  km  and  4>  is  the  latitude 

in  magnetic  coordinates.  In  this  expression  we  will  take  h^  to 

be  70  km  and  use  (34)  to  determine  sin  $.  The  angles  i],'^  and  are 

determined  from  the  relationships  i 


cos  tp . = sin 

e 

cos  cos 

I 

+ cos  0 

sin  I 

(37a) 

cos  ip  = sin 
r 

e 

cos  B cos 
r 

I 

- cos  0 

sin  I 

(37b) 

where  is  the  bearing  of  the  receiver  as  seen  from  the  transmitter 
in  magnetic  coordinates  and  I,  the  angle  of  inclination  of  the  earth's 
magnetic  field,  is  given  by 


10 


r 


tan  I = -2  tan  <I).  (38) 

For  multiple-hop  paths,  the  expression  (35)  must  be  multiplied 
by  the  number  of  hops. 

B. 2 Antenna  Gain 

No  attempt  has  been  made  in  this  model  to  provide  a variety 

of  antenna  patterns.  While  any  desired  antenna  gain  can  be  easily 

included  in  the  program,  we  shall  use  the  gain  of  a short  vertical 

dipole  for  both  transmitter  and  receiver.  If  G and  G denote 

T R 

the  transmitter  and  receiver  power  gains  relative  to  isotropic 
antennas  we  have 


G^(0)  = (0) 


3 . 2 
2 sin 


0 . 


Thus  for  the  overall  attenuation  in  decibels  we  have 


(39) 


B. 3 Spreading  Loss 

The  decrease  in  signal  strength  attributable  to  the  length  of 
the  transmission  path  is  given  by 

- 20  log^o  (41) 

where  is  the  path  loss  in  decibels,  D is  the  total  path  length 
in  km,  taking  into  account  any  multiple  hops,  and  f is  the  fre- 
quency in  mHz. 

The  total  attenuation  for  a particular  ionosphere  return  is 
obtained  by  adding  (35),  (40)  and  (41),  after  making  any  necessary 

adjustments  for  multiple  hops.  Denoting  the  total  attenuation  for 
the  j-th  return  by  we  have 

S * ‘•ant  * ‘-d-  l‘'2) 


‘1 


! 


li 
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The  rms  gain  is  thus  given  by 


A,  = . 

D 


The  remaining  quantities  needed  to  characterize  the  random 
gains,  (t) , are  the  doppler  shifts  and  the  doppler  spreads.  In 
this  report  we  are  not  concerned  with  the  relatively  large  doppler 
values  that  occur  during  sunrise  and  sunset,  but  with  the  smaller 
values  that  are  typical  of  the  ionosphere  for  the  majority  of  a 
typical  twenty-four  hour  period.  What  is  desired  is  the  shift  and 
spread  associated  with  each  of  the  propagation  modes.  Unfortunately, 
measurements  that  have  been  made  concerning  doppler  values  usually 
have  not  been  concerned  with  the  shift  and  spread  of  the  individual 
propagation  modes;  instead  they  have  recorded  composite  doppler 
values  involving  many  modes  [18].  Moreover,  most  theoretical 
studies  have  not  been  concerned  with  communication  system  per- 
formance and  are  of  marginal  value  for  present  purposes  [19-22]  . 

What  is  required  for  our  model  is  an  expression  relating  the 
doppler  shifts  and  spreads  to  frequency  for  each  ionosphere  return. 
This  expression  will  depend  on  the  electron  density  profile  and 
its  time  rate  of  change  as  well  as  on  various  other  ionosphere 
parameters  [23]  . VJe  assurt^e  that  the  doppler  values  can  be  char- 
acterized adequately  by  thp  following  power  laws 

V = A f°‘  (44) 


a = B f 


where  v and  2o  are  the  doppler  shift  and  doppler  spread;  f is  the 
frequency  in  mHz;  and  A,  B,  a,  and  6 are  constants  that  characterize 
the  region  from  which  the  wave  is  reflected  and  the  magnetoionic 
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component  involved.  In  general  we  need  an  expression  of  the  above 
type  for  the  low  and  high  ray  returns  from  each  region  and  for  each 
magnetoionic  component. 

It  is  convenient  to  express  the  constants  A and  B appearing  in 
(44)  and  (45)  in  terms  of  a reference  value  measured  at  a known 
frequency.  Thus  the  above  expressions  can  be  rewritten  as 


V = v^(f/fj.)^  (46) 

a = a^(f/fj.)^  (47) 

where  and  2a^  are  the  doppler  shift  and  spread  measured  at  the 
reference  frequency  f^..  For  multiple  hops  we  assume  that  the  doppler 
effects  combine  incoherently  so  that 


and 


V = (No.  of  Hops)  • (f/f  j.)^ 

U 6 

a = (No.  of  Hops)  (f/f  j.) 


(48) 

(49) 


Given  the  appropriate  reference  values  and  exponents,  the  above 
expressions  provide  the  doppler  shift  and  spread  for  each  ionosphere 
return.  However,  since  the  reference  values  and  exponents  usually 
are  not  known,  we  must  make  some  assumptions  to  obtain  "default" 
values  for  use  in  the  absence  of  more  accurate  data.  It  will  be 
assumed  that  the  doppler  shift  and  spread  are  the  same  for  the  high 
ray  and  low  ray  and  for  both  magnetoionic  components  of  the  E and 
F region  returns.  The  following  table  gives  the  "default"  reference 
values  and  exponents  to  be  used  for  the  E and  F regions.  These 
values  are  based  upon  measurements  made  by  Watterson,  et.  al.  [1]. 
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Region 

Vj.(Hz) 

2o^  (Hz) 

f j.  (mHz) 

a 

E 

0.01 

0.02 

9.30 

1.0 

1.0 

F 

0.01 

0.15 

9.30 

1.0 

1.0 

Table  1.  Default  Values  for  Doppler  Shift  and 
Doppler  Spread 

D.  Computer  Program 

The  Fortran  program  used  to  calculate  the  foregoing  channel 
parameters  is  presented  in  Appendix  B.  The  basic  output  of  this 
program  is  a listing  of  the  input  parameters  and  a tabular  presen- 
tation of  the  computed  channel  parameters.  A typical  output  is 
shown  in  Table  2.  The  indexing  convention  used  for  the  returns  is 
presented  in  Table  3 and  a glossary  of  the  terms  used  in  the  out- 
put is  presented  in  Table  4. 

A modified  version  of  this  output  can  be  produced  in  which  all 
returns  whose  attenuation  exceeds  some  threshold  value  above  the 
minimum  attenuation  are  omitted.  Those  modes  for  which  there  is  no 
return  are  likewise  omitted.  Using  this  modified  output  one 
can  identify  those  returns  that  are  of  primary  interest  for 
communication  purposes.  An  example  is  shown  in  Table  5 using 
the  same  inputs  that  produced  Table  2 and  an  attenuation  threshold 
of  40  dB. 

A synthetic  ionogram  also  has  been  produced  using  the  program. 
The  resulting  curves  are  presented  in  Figure  5 and  include  returns 
for  one,  two  and  three  hops  together  with  the  corresponding  attenu- 
ations . 
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Table  3 


Index  Convention  for  Ionosphere  Returns 


Index 

Number 
of  Hops 

Region 

Ordinary  or 
Extraordinary 

Low  Ray  or 
High  Ray 

1 

1 

E 

0 

LOW 

2 

1 

E 

0 

HIGH 

3 

1 

F 

0 

LOW 

4 

1 

V 

0 

HIGH 

5 

1 

E 

X 

LOW 

6 

1 

E 

X 

HIGH 

7 

1 

F 

X 

LOW 

8 

1 

F 

X 

HIGH 

9 

2 

E 

0 

LOW 

10 

2 

E 

0 

HIGH 

11 

• 

2 

• 

F 

• 

0 

• 

LOW 

• 

• 

44 

• 

6 

• 

F 

• 

0 

• 

HIGH 

45 

6 

E 

X 

LOW 

46 

6 

E 

X 

HIGH 

47 

6 

F 

X 

LOW 

48 

6 

F 

X 

HIGH 
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Table  4 . Glossary  of  Terms  Used 

for  Computed  Channel  Parameters 


Mode 

Identifies  particular  ionosphere 

return  (e.g.,  2F0  LOW  is  the  two- 

hop,  F region,  ordinary  wave,  low 
ray  return) 

Solution  Indicator 

This  equals  1 if  there  is  a return 

for  the  particular  mode  and  equals 

0 if  there  is  no  return 

Ray  Angle 

1 

0 

1 

Total  Path  Length 

D . 

1 

Carrier  Delay 

T(fc) 

Carrier  Phase 

The  fractional  part  of  D.^°^ 

Signal  Delay 

D.(l) 

D 

Amplitude  Distortion 

D.^2) 

: 

Phase  Distortion 

d.(3> 

D 

Attenuation 

L . 

I 

Doppler  Shift 

V . 

D 

Doppler  Spread 

2a  . 

: 
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IV.  SIGNAL  ANALYSIS 


This  section  is  concerned  with  the  signal  analysis  part  of  the 
simulator  and  the  approximations  used  in  its  implementation. 

We  will  assume  that  the  transmitted  signal,  s(t),  is  bandlimited 
with  (two  sided)  bandwidth  BW.  The  received  signal,  r(t),  is  not 
strictly  bandlimited  because  it  is  multiplied  by  the  random  gains, 
which  are  not  bandlimited.  However,  since  the  doppler  spreads  are 
small  compared  with  the  signal  bandwidth,  the  received  signal  can 
also  be  assumed  to  be  bandlimited  but  with  a slightly  larger  band- 
width than  the  transmitted  signal.  The  error  introduced  by  this 
assumption  is  negligible  as  long  as  the  assumed  bandwidth  of  the 
received  signal  is  large  compared  with  the  sum  of  the  transmitted 
signal  bandwidth  and  the  largest  doppler  spread  present.  In  the 
following  we  will  assume  that  all  signals  have  a bandwidth  equal 
to  twice  that  of  the  transmitted  signal.  This  is  sufficient  to 
reduce  aliasing  errors  to  a negligible  value.  Denoting  the 

sampling  interval  associated  with  this  larger  bandwidth  by  T we  have 

o 


T = ^ 
o 2*BW 


(50) 


As  a consequence  of  the  cardinal  series  representation  for  band- 
limited  functions,  we  require  a knowledge  of  the  signals  involved 
only  at  the  sample  times  nT^,  n = 0,  ±1,  •••.  This  greatly  simpli- 
fies the  computational  and  storage  requirements.  From  (21)  the 
received,  sampled  signal  appears  as 


r (nT^) 


N 

= A s(nT  -N  T ) + S exp[-2TTiD . 

Q O O O . , 1 

:=i 


•G . (nT  ) s . (nT  -N  T ) 
1 o : o j o 


n=0 ,!,••• 


(51) 


where  the  various  delays  have  been  replaced  by  the  nearest  integer 
multiple  of  T^ . That  is,  Ng  is  defined  as  the  integer  satisfying 


29 


(52) 


D - 4t  <NT  <D  + iT 
9 ^ o — g o g 2o 


and  the  N.'s  are  defined  as 
D 


the  integers  satisfying 


D . 


(1) 


1 T < M . T D . 

2 o - 3 o 3 


j = lr 


•N. 


(53) 


It  remains  to  determine  the  sampled,  distorted  signals 

and  the  sampled  random  path  gains  G.(nT  ). 

3 o 


s . (nT  ) 
3 o 


A.  Delay  Distortion 

The  transmitted  signal  has  the  representation 


s ( t) 


sin  2nBW(t-nT  ) 

s (nT  ) — 2_ 

o 27iBW(t-nT  ) 


(54) 


and  the  impulse  response  h^ (t) , which  we  also  assume  to  be  band- 
limited,  has  the  representation 


hj  (t) 


oo 

h . (nT  ) 
3 o 

n=-oo 


sin  27TBW(t-nT  ) 
Q_ 

2TTBW(t-nT  ) 
o 


j=l,---N.  (55) 


V'?e  obtain  the 
hj (t) . But  Sj 


distorted  waveform,  s^  (t) , by  convolving  s(t)  and 
(t)  also  has  an  expansion  of  the  above  type  and 


^ -s  2’nBW(t-nT  ) 

s.(t)  = 7 s.(nT  ) ^ 

^ 3 o 2TTBW(t-nT  ) 

n=-oo  o 


3 


, j = l , ’ • -N.  (56) 


We  are  interested  here  in  the  sample  values  s . (nT  ) and  those  are 

3 o 

given  by  the  convolution 


s . (nT  ) 
3 o 


s(nT^-mT^)  h^  (mT^)  , j = l,---,N. 

m=-oo 


(57) 
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J 


where 


BW 


h.(mT^)  = T^h.(mT^) 


= J"  [2TTifinTQ  ] df . 

-BW 

j=l , • • ’N 


(58) 


The  transfer  function  IK  (f ) is  given  by  (11),  and  it  will  be 
approximated  with  the  first  two  terms  of  the  series  appearing  in 
the  exponent,  that  is 


Hj  (f ) ::  exp[-27rif  - 2Tiif  ^ ^ ^ ^ , j-l,**-N. 

After  the  change  of  variable  n = f/BW,  we  obtain  from  (58) 


hj  (mT^)  ~ j J'  exp[-iAj  • expIimTTnJd  n 

-1 

j=l, • • -N 


(59) 


(60) 


with 


A . = 2tt  (BW)  ^ D . 


Aj  = 2tt  (BW)  ^ 


(61) 

(62) 


Closed  form  expressions  for  (60)  are  not  known  and  further 
approximations  are  required  for  the  numerical  determination  of 
hj (mT^) . Let  I(x)  represent  the  general  integral  defining  the 
impulse  responses , that  is 


I(x)  - \ J"  exp[-iA^^^  n ^-iA  n^]*  exp[ixn]dr 
-1 


(63) 


(3)  3 

Expanding  exp[-iA  n ] in  a power  series  we  have  the  representation 
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I 


I 


IM 

p=0 


(-iA 


(3) 


P 

) 


P! 


I 


P 


(x) 


where  we  have  defined 


1 

lp(x)  = i J"  exp [-iA n^] ‘exp [ixn] dn 
-1 


(64) 


(65) 


The  integral  Iq(x)  is  readily  expressed  in  terms  of  Fresnel  inte- 
grals [24] . For  the  higher  order  integrals  the  following  recurrence 
relation  has  been  derived 


Ip(x)  = ^2)  exp[-iA^^^]  |exp[ix]  + (-1)^  exp[-ix]| 


^ ~ > I T (x)  -i  I 


(66) 


2A 


(2)  "p-1 


2A 


(2)  ^p-2 


where  we  have  adopted  the  convention  that  l_^(x)=0. 

The  impulse  responses  for  the  ionosphere  returns  can  now  be 
evaluated  with  relative  ease.  Starting  with  Iq(x)  and  using  (66) 
we  determine  as  many  of  the  terms  ip(x)  as  are  needed  to  adequately 
approximate  the  infinite  sum  (64).  The  integral  Iq(x)  is  evaluated 
using  the  rational  approximation  presented  in  [25]  . Care  must  be 
exercised  in  implementing  this  procedure,  however.  The  error  in- 
curred by  using  the  rational  approximation  for  I (x)  is  approximately 

-9  ° 

1.0  X 10  . For  large  arguments  (x>5)  nearly  equal  functions  are 

subtracted  in  (66)  and  this  reduces  the  number  of  significant 

figures  in  the  resulting  expression.  For  the  higher  order  terms 

this  reduction  in  the  number  of  significant  figures  is  important 

and  the  computed  values  may  become  inaccurate.  It  is  evident  from 

(64)  that  this  effect  will  be  more  noticeable  when  A^^^  is  large. 

Experience  has  shown  that  values  of  x as  large  as  20-25  can  be 

tolerated  in  most  applications  before  the  errors  become  excessive. 


I 

i 


I 
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B. 


Random  Gains 


The  random  gains,  (t) , in  the  ionosphere  returns  are  complex 


valued,  zero  mean,  Gaussian  processes  whose  correlation  function  and 
power  spectral  density  are  given  by  (cf.  Eqns.  19  and  20) 


(t)  = A.^  exp  [-2TT^a . + 2iTiv.T]  , j=l, 

j ^ : 


■N 


(67) 


and 


A. 

(f)  = — ^ 


/2TTa  . 
J 


exp 


r ^ 

I 2“/ 


2 


f j=lf 


•N. 


(68) 


We  wish  to  generate  a random  sequence  with  the  above  correlation 
function  and  we  wish  to  do  this  directly  in  the  time  domain  and  not 
rely  on  Fourier  transform  techniques.  Because  the  doppler  spreads. 


Oj,  are  small  relative  to  the  signal  bandwidths,  the  time  domain 


approach  requires  a large  memory.  To  avoid  these  large  memory 
requirements  we  will  replace  the  original  process  with  an  approxi- 
mation that  can  be  generated  more  easily. 


The  approximate  random  gain  will  be  denoted  by  G^  (t)  and  it 


will  be  chosen  to  have  the  power  spectral  density 

2 


8A.' 


p;^  (f)  = - 


» j=l» 


/N 


(69) 


and  corresponding  correlation  function 

(t)  = Aj^  1^1  + 2TTpaj|T|  + j ( 2iTpaj  I T I )^] 


•exp 


-2Tr 


POj 


+ 2 


TTiVjT  , j=l,*--. 


(70) 


N, 


The  parameter  p is  chosen  so  that  (70)  agrees  with  (67)  over  a 
reasonable  range  of  arguments.  The  value  p =2.1  has  been  deter- 
mined to  be  a good  choice  and  this  is  the  value  that  will  be  used 


>1 
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throughout  the  following.  Normalized  plots  of  C„  (x)  and  (i) 
are  compared  in  Figure  6 . ^ ^ 

It  is  obvious  from  (69)  that  the  approximate  process  is 

not  bandlimited.  However,  the  error  introduced  by  assuming  that 
the  process  has  bandwidth  BW  will  be  negligible  since  (o./BW)<<l, 
Accordingly  we  will  treat  (t)  as  if  it  is  bandlimited,  in  which 
case  we  need  only  evaluate  it  at  the  sample  times  nT^.  Thus,  the 
original  problem  has  been  reduced  to  the  generation  of  a complex 
valued,  zero  mean,  Gaussian  sequence  with  correlation  function 


Cv  (nT  ) 

Gj  ° 


A.^l^l  + 2Trpa.T^|n|  + i ( 27Tpa  . | n | ) ^ J 

•expj^-2TrpajTQ|  n I + 27riVjTQnl  . 

n=0,±l,... 
j=l , • • ‘N 


(71) 


The  technique  used  for  the  computer  generation  of  the  approximate 
process  is  discussed  in  Appendix  C.  A plot  of  the  normalized, 
empirical  correlation  function  of  the  computer-generated  process 
is  also  displayed  in  Figure  6.  Agreement  between  the  computer- 
generated correlation  function  and  the  approximate  correlation 
function  (71)  is  excellent  for  the  range  of  arguments  considered. 

C . Computer  Program 

A listing  of  the  Fortran  program  used  to  implement  the  s^ignal 
analysis  part  of  the  simulator  is  presented  in  Appendix  D. 

Channel  parameter  information  is  read  into  the  program  via  data  cards 
that  have  been  produced  by  the  first  part  of  the  simulator.  The 
sampled  input  values  are  read  from  a tape,  the  program  calculates 
the  sampled  output  values  according  to  Eqn.  (51)  and  writes  the 
output  on  another  tape.  As  an  alternate,  printed  output  also  can 
be  obtained,  but,  because  the  quantity  of  output  data  is  large,  the 
usefulness  of  this  mode  of  presentation  is  severely  limited. 


V.  CONCLUDING  REMARKS 


RECOMMENDATIONS 


The  simulator  described  in  the  foregoing  sections  represents 
a significant  improvement  over  earlier  simulators.  It  has  a 
wideband  capability  not  previously  available  and  its  inputs  are  the 
basic  physical  quantities  that  characterize  the  ionosphere  and  the 
transmitter-receiver  environment.  Explicit  values  for  the  channel 
parameters  are  not  required  as  these  are  calculated  internally 
using  the  model  of  the  ionosphere.  Comparison  of  the  computer- 
generated ionogram.  Figure  5,  and  ionograms  actually  observed  in 
practice  [27] , suggests  that  our  mathematical  model  is  adequate 
for  most  communication  purposes.  However,  a precise  comparison  of 
the  simulator  and  an  actual  channel  is  difficult  because  of  the 
variability  and  lack  of  accuracy  of  much  of  the  available  iono- 
spheric data. 

Two  aspects  of  the  present  simulator  deserve  further  study. 

The  first  is  concerned  with  the  absorption  mechanism,  while  the 
second  is  concerned  with  the  determination  of  appropriate  input 
parameters.  The  absorption  predicted  for  the  F region  high  ray 
return  is  noticeably  smaller  than  the  absorption  usually  observed 
in  practice.  Since  both  the  delay  and  delay  distortion  of  the 
high  ray  is  large,  this  can  have  an  adverse  affect  on  signal  trans- 
mission. In  actual  practice  these  returns  suffer  a large  attenu- 
ation, and  their  influence  is  relatively  small.  However,  our 
model  does  not  provide  enough  absorption  to  reduce  the  high  ray 
signal  strengths  to  negligible  levels. 

The  remedy  may  be  simply  a matter  of  choosing  different  con- 
stants in  the  nondeviative  absorption  expression,  although  prelimin- 
ary results  do  not  fully  support  this  approach.  More  likely  it  will 
be  necessary  to  introduce  a deviative  absorption  term  in  the 
equation  for  total  loss.  The  inclusion  of  such  a term  would  also 
reduce  the  signal  strengths  for  nighttime  transmission,  and  this 
too  would  make  the  total  loss  agree  more  closely  with  observations. 

The  second  area  for  further  study  is  concerned  with  determining 
suitable  input  parameters  so  that  the  performance  predicted  by  the 


mathematical  model  agrees  as  well  as  possible  with  the  actual 
observed  performance.  This  involves  matching  actual  and  computer- 
generated observables  such  as  doppler  shift,  doppler  spread,  attenu- 
ation and  delay  vs.  frequency  plots  (ionograms)  to  determine  the 
input  pareimeters.  Besides  the  obvious  advantage  of  providing  further 
justification  of  the  assumptions  used  in  the  model,  this  parameter 
determination  would  provide  the  user  with  a set  of  inputs  that 
correspond  to  a variety  of  actual  environments. 

The  primary  difficulty  encountered  in  the  determination  of 
the  input  parameters  is  that  an  explicit  functional  relationship 
does  not  exist  between  the  input  parameters  and  the  computed  channel 
parameters;  the  channel  parameters  are  determined  by  numerically 
solving  a transcendental  equation.  As  a consequence,  techniques 
must  be  developed  that  are  tailored  to  the  ionosphere  model.  One 
approach  is  to  select  the  most  significant  observable  features  and 
to  determine  the  relationship  between  these  features  and  the  input 
parameters.  For  example,  the  nearly  horizontal  portions  of  an 
observed  ionogram  would  be  affected  most  by  the  heights  of  the  regions 
and  to  a lesser  degree  by  the  thickness  of  the  regions  and  their 
peak  electron  densities.  On  the  other  hand,  the  maximum  operating 
frequency  is  likely  to  be  more  affected  by  the  peak  electron  densi- 
ties than  the  heights  or  thicknesses  of  the  regions.  A detailed 
study  of  these  and  other  pertinent  relationships  would  provide  the 
user  with  the  ability  to  determine  input  parameters  corresponding 
to  a variety  of  observed  operating  conditions.  These  same  condi- 
tions could  then  be  reproduced  by  the  channel  simulator  as  required 
for  system  studies. 
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APPENDIX  A 

DELAY  AND  DELAY  DERIVATIVES 


In  this  appendix  we  present  the  equations  for  the  delay  and 
the  first  three  delay  derivatives.  Recalling  the  definition  of 
A(x)  we  have 

o *5 


A(x)  = 

where  x = cos  0 . VJe  must  solve  the  equation 


(A.l) 


R = A(x) 


(A. 2) 


with  R equal  to  the  distance  between  transmitter  and  receiver.  The 
solution  then  yields  the  delay  via  the  equation 


T = 
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with  c the  velocity  of  light.  The  delay  derivatives  are  given  by 


(A. 3) 


dr 

df 


XT  dx 


(A. 4) 


2 2 
d T l+2x 


df 


2 

X T 


(i) 


XT  d X 


(1-x  ) df' 


(A. 5) 


and 


d^T 


3 (l+x^+2x‘^) 


df 


4 2 

X T 


(IJ) 


^ ^ 3(l+2x^) 


2 

X T 


df' 


XT  d X 


(A. 6) 


(l-x*-)  df 


A-1 


■ST 


I ■.  IWWMmM 


Thus  to  determine  the  delay  and  delay  derivatives  we  must  solve 
(A. 2)  for  X and  then  evaluate  the  first  three  derivatives  of  x 
with  respect  to  f for  use  in  (A. 4)  - (A. 6).  In  the  follov/ing 
these  functions  are  listed  for  the  two  electron  density  profiles. 
The  group  height  G(x)  is  obtained  by  evaluating  the  integral  (24). 

1.  Nighttime 

For  the  nighttime  electron  density  the  integral  for  the  group 
height  is  evaluated  and  used  in  (A.l)  to  yield 
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For  fj,x^l  the  derivatives  of  x with  respect  to  frequency  are: 
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2 . Daytime 

The  daytime  electron  density  profile  differs  from  the  night- 
time profile  only  for  heights  greater  than  h . Thus  the  daytime 
equations  for  f„<l  are  the  same  as  the  nighttime  equations. 

In  this  case 
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where 


r = 


_pE 


pF 


For  fpX£l  the  derivatives  of  x with  respect  to  frequency  are  given 


by  (A. 8)  - (A. 10).  For  fpX>l  we  have 
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PROGRAM  FOR  CALCULATING  CHANNEL  PARAMETERS 
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APPENDIX  C 

COMPUTER  GENERATION  OF  RANDOM  GAINS 


For  each  ionosphere  return,  it  is  desired  to  generate  a 
stationary,  zero-mean  Gaussian  sequence  (nT^) , n = 0,1,---, 
having  the  correlation  function 


eg  (nT^) 


= E (nTg  + mTo)  (mT^) 

= (27TpajTQ  |n  I ) + j (27TpajTQ|n|)^J 

•exp  |^-2iTpa  jT^  I n I + 2iriVjT^nj  (C 


(C.l) 


n = 0 , 1 , • • • 
j = 1, • • • ,N  . 

We  will  generate  this  sequence  by  passing  a white  noise  process 
through  an  appropriately  chosen  discrete-time  linear  filter  [26] . 
Denoting  the  filter  input  by  {Cj(n)}  and  the  impulse  response  by 
{a^ (n) } we  have 


(nTo)  = exp[2TTiVjnTQ]^^  5j(n-m)aj 


(C.2) 


j=l, • • ’N 


where 


a.(m)  = 

^ 0 , 


AjdirpOjTo)'^  (2iTpa  jmTol^exp  [^-2TTpajmT  J , m>0 


(C.3) 


j=l, • • -N. 


C-  1 


The  random  variables  are  independent,  identically-distributed, 

complex  Gaussian  variables.  They  may  be  expressed  as 


Cj(n)  = + i^PNn)  , n=0,l,...  (C.4) 

where  ^f^Nn)  and  C-^Nn)  are,  respectively,  the  real  and  imaginary 
J 3 n ) ( 2 ) 

parts  of  Cj  (n)  . For  each  n the  variables  ? j (n)  and  f; ^ (n)  are 

independent,  zero-mean,  unit  variance,  real  Gaussian  variables,  and 

both  components  of  (n)  are  independent  of  each  other  for  all  values 

of  n . 

To  generate  G^ (nT^)  directly  using  the  convolution  (C.2), 
we  would  have  to  approximate  the  infinite  summation  with  a sufficiently 
accurate  truncated  summation.  Because  OjT^  is  very  small,  the  length 
of  this  truncated  summation  would  be  large  and  we  would  expect  error 
accumulation  and  an  attendant  loss  of  accuracy.  However,  our  par- 
ticular choice  of  an  approximating  correlation  function  (C.l)  allows 
us  to  generate  the  random  sequence  G^ (nT^)  iteratively.  This  was, 
of  course,  the  primary  consideration  in  selecting  this  particular 
approximation.  The  iterative  scheme  is  given  below 


G.  (nT  + 3T  )=  K.  exp[-2TTpa.T  + 2Triv . (n+3 ) T ] ^ . (n+2) 

j o o J JO  3 o;] 

+ K.  exp  [-4'n-pa  .T  + 2TTiv.(n+3)T  ] f.(n+l) 

j JO  J O ] 

+ exp  [-6iTpa jT^  + eiriVjT^]  (nT^)  (C.5) 

- 3 exp[-47rpa.T  + 4Triv.T  ] (nT  + T ) 

*"  Do  D03  0 o 

+ 3 exp [-2Trpa jT^  + 2iTiVjT^l  ^ ^"^o^ 

n>0 

j=l , • • ’N 


1 


C-2 


where 


L 4 


% 

K.  = A.  (2TTpa  T ) , j = l,---N 

^ /3  3 J 

and,  as  discussed  earlier,  we  will  choose  p = 2.1.  To  start  this 
iteration  we  need  0^(0),  (T^)  and  ^j(2T^)  given  by 

(0)  = AjWj  (1) 


(C.6) 


S . (T  ) = exp[2iTiv  .T  ] • Y ■ (1)  W (2)]  + X . (1)  S . (0)  1 
jo  joD  j j j* 

^j(2T^)  = exp  [4iriVjT^]  • j Yj  (2)  Wj  (3)  + exp  [-2ttpo  jT^]  (1 ) 

+ X . (1)  [1-A  . (2)]  (T  ) + X.(2)  &.(0)N 
3 3 3 0 3 3 ' 


(C.8) 


(C.9) 


j=l, ’ ' 'N. 


Here  (1) , (2) , (3) , and  (1)  are  independent,  complex, 

Gaussian  variables  v;hose  real  and  imaginary  parts  are  independent 
and  have  zero  mean  and  unit  variance.  The  other  terms  appearing 
in  (C.7)-(C.9)  are 

1 2 

X.(l)  = [1+(4ttpo.T  ) + ^ (4iTpa.T  ) ] ‘exp  [-4TTpa  .T  ] 

3 3o3  30  30 

X.(2)  = [1+(8ttpo.T  ) + i (Strpo.T  ) ^]  ‘exp  [-Sirpo  .T  ] 

3 J o 3 JO  JO 

Yj(l)  = A.  [l-X^j(l)]^ 

Yj(2)  = Aj  l^^l-X^  j (1)^  ^1-X^  j (2)^  - j ^2TTpOjT^^^  exp  [-4ttpo  ^T^lj  ' 


(C.IO) 


(C.ll) 


(C.12) 


(C.  13) 


The  coefficients  A^  are  determined  by  the  attenuation  of  each  return 
via  Eqn.  (43) . 


It  remains  to  generate  the  sequence  of  independent,  Gaussian 
variables  {^^(n)}.  We  begin  by  using  the  congruence  method  to 
generate  a sequence  of  uniform  random  variables.  This  method 
generates  a sequence  of  random  integers  via  the  iterative  scheme 
[24] 


Y = a'Y  [mod  M]  . 
n+1  n 


(C.14) 


5 *5  1 

For  the  quantities  a and  M,  we  have  found  a = 7 and  M = 2 to  be 
effective.  Multiplication  is  done  using  double  precision  arithmetic 
but  the  final  result  is  expressed  as  a single  precision  variable. 
Upon  division  of  {Y  } by  2^^  a sequence  of  independent,  uniform 
[0,1]  variables  is  produced.  This  sequence  is  then  transformed 
into  a sequence  of  independent,  complex,  Gaussian  variables. 

If  u^  and  U2  are  real,  independent,  uniform  variables,  the  real 
and  imaginary  parts  of  the  variable  E,  defined  by 


C = (-2  1nUj^)^  exp[2iTiu2] 


(C.15) 


are  independent  and  Gaussian  with  zero  mean  and  unit  variance  [24]. 
Thus,  we  obtain  the  desired  sequence  of  complex  Gaussian  variables 
for  use  as  the  filter  input  in  Eqn.  (C.5). 
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